HEALPix logo

UXarray for Basic HEALPix Statistics & Visualization

In this section, you’ll learn:

  • Utilizing intake to open a HEALPix data catalog

  • Using the uxarray package to look at basic statistics over HEALPix data

  • Using UXarray plotting functionality on HEALPix data

Prerequisites

Concepts

Importance

Notes

UXarray

Necessary

HEALPix overview

Necessary

Time to learn: 30 minutes

import uxarray as ux
import cartopy.crs as ccrs

Open data catalog

Tip

If you think you first need to learn about Intake, Pythia’s Intake Cookbook is a great resource to do so.:::

Tip

Also, if you want to get insights into the particular HEALPix data catalog we will use throughout this section, please be sure to check the previous section first, easy.gems for HEALPix Analysis & Visualization:::

import intake

cat = intake.open_catalog("https://data.nextgems-h2020.eu/online.yaml")
model_run = cat.ICON.ngc3028

ds_coarsest = model_run().to_dask()
ds_fine = model_run(zoom=7).to_dask()

/home/runner/miniconda3/envs/healpix-cookbook-dev/lib/python3.10/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

/home/runner/miniconda3/envs/healpix-cookbook-dev/lib/python3.10/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

uxds_coarsest = ux.UxDataset.from_healpix(ds_coarsest)

uxds_fine = ux.UxDataset.from_healpix(ds_fine)
uxds_fine

<xarray.UxDataset> Size: 2TB
Dimensions:                              (time: 1350, depth_half: 129,
                                          n_face: 196608, level_full: 90,
                                          crs: 1, depth_full: 128,
                                          soil_depth_water_level: 5,
                                          level_half: 91,
                                          soil_depth_energy_level: 5)
Coordinates:
  * crs                                  (crs) float32 4B nan
  * depth_full                           (depth_full) float32 512B 1.0 ... 5....
  * depth_half                           (depth_half) float32 516B 0.0 ... 6....
  * level_full                           (level_full) int32 360B 1 2 3 ... 89 90
  * level_half                           (level_half) int32 364B 1 2 3 ... 90 91
  * soil_depth_energy_level              (soil_depth_energy_level) float32 20B ...
  * soil_depth_water_level               (soil_depth_water_level) float32 20B ...
  * time                                 (time) datetime64[ns] 11kB 2020-01-2...
Dimensions without coordinates: n_face
Data variables: (12/88)
    a_tracer_v_to                        (time, depth_half, n_face) float32 137GB ...
    atmos_fluxes_frshflux_evaporation    (time, n_face) float32 1GB ...
    atmos_fluxes_frshflux_precipitation  (time, n_face) float32 1GB ...
    atmos_fluxes_frshflux_runoff         (time, n_face) float32 1GB ...
    atmos_fluxes_frshflux_snowfall       (time, n_face) float32 1GB ...
    atmos_fluxes_heatflux_latent         (time, n_face) float32 1GB ...
    ...                                   ...
    va                                   (time, level_full, n_face) float32 96GB ...
    vas                                  (time, n_face) float32 1GB ...
    w                                    (time, depth_half, n_face) float32 137GB ...
    wa_phy                               (time, level_half, n_face) float32 97GB ...
    wind_speed_10m                       (time, n_face) float32 1GB ...
    zos                                  (time, n_face) float32 1GB ...

HEALPix basic stats using UXarray

import matplotlib.pylab as plt

boulder_lon = -105.2747
boulder_lat = 40.0190

Mesh face containing Boulder’s coords

%%time
boulder_face = uxds_fine.uxgrid.get_faces_containing_point(point_lonlat=[boulder_lon, boulder_lat])

CPU times: user 8.46 s, sys: 96 ms, total: 8.55 s
Wall time: 8.42 s

Global and Boulder’s temperature averages

In order to get a line plot out of our UXarray.UxDataset objects’ 1-dimensional temperature variables, we will convert them to xarray and call the default plot function because UXarray’s default plot functions are all dedicated to grid-topology based visualizations:

%%time
uxds_fine.tas.isel(n_face=boulder_face).to_xarray().plot(label="Boulder")
uxds_coarsest.tas.mean("n_face").to_xarray().plot(label="Global mean")

plt.legend()

CPU times: user 3.56 s, sys: 1.09 s, total: 4.65 s
Wall time: 6.9 s

<matplotlib.legend.Legend at 0x7eff02f609d0>
../_images/450dbbbf1e6fbb0967b8d6ddd460b3384fb03e1e9db6489be091db043052984f.png

Data plotting with UXarray

UXarray provides several built-in plotting functions to visualize unstructured grids, which can also be applied to HEALPix grids in the same interface:

Let us first look into interactive plots with bokeh backend (i.e. UXarray’s plotting functions have a backend parameter that defaults to “bokeh”, and it can also accept “matplotlib”)

Global plots

%%time

projection = ccrs.Robinson(central_longitude=-135.5808361)

uxds_fine["tas"].isel(time=0).plot(
    projection=projection,
    cmap="inferno",
    features=["borders", "coastline"],
    title="Global temperature",
    width=700,
)

CPU times: user 7.36 s, sys: 28.4 ms, total: 7.39 s
Wall time: 7.96 s

WARNING:param.GeoOverlayPlot00477: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

%%time

uxds_fine["tas"].isel(time=0).plot(
    backend="matplotlib",
    projection=projection,
    cmap="inferno",
    features=["borders", "coastline"],
    title="Global temperature",
    width=1100,
)

CPU times: user 485 ms, sys: 3.08 ms, total: 488 ms
Wall time: 1.71 s

Regional subsets (Not only for plotting but also for analysis)

When a region on the globe is of interest, UXarray provides subsetting functions, which return new regional grids that can then be used in the same way a global grid is plotted.

Let us look into USA map using the Boulder, CO, USA coords we had used before for simplicity:

Subsetting uxds_fine into a new UxDataset using a “bounding box” around Boulder, CO first:

%%time

lon_bounds = (boulder_lon-20, boulder_lon+40)
lat_bounds = (boulder_lat-20, boulder_lat+12)

uxds_fine_subset = uxds_fine["tas"].isel(time=0).subset.bounding_box(
    lon_bounds, 
    lat_bounds
    )

/home/runner/miniconda3/envs/healpix-cookbook-dev/lib/python3.10/site-packages/uxarray/grid/grid.py:1432: RuntimeWarning: Necessary functions for computing the bounds of each face are not yet compiled with Numba. This initial execution will be significantly longer.
  warn(

CPU times: user 1min 5s, sys: 107 ms, total: 1min 5s
Wall time: 1min 3s

If we quick check the global and regional subset’s average temperature at the first time-step, we can see the difference:

print("Regional subset's temperature average: ", uxds_fine_subset.mean("n_face").values, " K")
print("Global temperature average: ", uxds_fine.tas.mean("n_face").isel(time=0).values, " K")

Regional subset's temperature average:  278.97748  K

Global temperature average:  286.2222  K

Now, let us plot the regional subset UxDataset:

%%time 
projection = ccrs.Robinson(central_longitude=boulder_lon)

uxds_fine_subset.plot(
    projection=projection,
    cmap="inferno",
    features=["borders", "coastline"],
    title="Boulder temperature",
    width=1100,
)

CPU times: user 67.7 ms, sys: 0 ns, total: 67.7 ms
Wall time: 67.4 ms

Grid topology exploration

Exploring the grid topology may be needed sometimes, and UXarray provides functionality to do so, both numerically and visually. Each UxDataset or UxDataArray has their associated Grid object that has all the information such as spherical and cartesian coordinates, connectivities, dimensions, etc. about the topology this data belongs to. This Grid object can be explored as follows:

uxds_fine.uxgrid

<uxarray.Grid>
Original Grid Type: HEALPix
Grid Dimensions:
  * n_node: 196610
  * n_face: 196608
  * n_max_face_nodes: 4
Grid Coordinates (Spherical):
  * node_lon: (196610,)
  * node_lat: (196610,)
  * face_lon: (196608,)
  * face_lat: (196608,)
Grid Coordinates (Cartesian):
  * node_x: (196610,)
  * node_y: (196610,)
  * node_z: (196610,)
Grid Connectivity Variables:
  * face_node_connectivity: (196608, 4)
Grid Descriptor Variables:
  * n_nodes_per_face: (196608,)

There might be times that the user wants to open a standalone Grid object for a HEALPix grid (or any other unstructured grids supported by UXarray) without accessing the data yet, which can also be achieved as follows:

# Let us open a pretty coarse HEALPix grid because we will visualize it later
uxgrid = ux.Grid.from_healpix(zoom=3, pixels_only=False)
uxgrid

<uxarray.Grid>
Original Grid Type: HEALPix
Grid Dimensions:
  * n_node: 770
  * n_face: 768
  * n_max_face_nodes: 4
Grid Coordinates (Spherical):
  * node_lon: (770,)
  * node_lat: (770,)
  * face_lon: (768,)
  * face_lat: (768,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (768, 4)
Grid Descriptor Variables:

uxgrid.plot(
    periodic_elements="split",         
    projection=ccrs.Mollweide(),
    width=800,
    title="HEALPix, zoom=3")
On this page